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Abstract 

Suppose that we observe noisy linear measurements of an unknown signal that can be mod- 
eled as the sum of two component signals, each of which arises from a nonlinear sub-manifold of 
a high-dimensional ambient space. We introduce SPIN, a first-order projected gradient method 
to recover the signal components. Despite the nonconvex nature of the recovery problem and the 
possibility of underdetermined measurements, SPIN provably recovers the signal components, 
provided that the signal manifolds are incoherent and that the measurement operator satisfies 
a certain restricted isometry property. SPIN significantly extends the scope of current recovery 
models and algorithms for low-dimensional linear inverse problems and matches (or exceeds) 
the current state of the art in terms of performance. 

1 Introduction 

Estimation of an unknown signal from linear observations is a core problem in signal processing, 
statistics, and information theory. Particular energy has been invested in problem instances where 
the available information is limited and noisy and where the signals of interest possess a low- 
dimensional geometric structure. Indeed, focused efforts on certain instances of the linear inverse 
problem framework have spawned entire research subfields, encompassing both theoretical and 
algorithmic advances. Examples include signal separation and morphological component analysis [1, 
2]; sparse approximation and compressive sensing [3-5]; affine rank minimization [6]; and robust 
principal component analysis [7, 8] . 

In this work, we study a very general version of the linear inverse problem. Suppose that the 
signal of interest x* can be written as the sum of two constituent signals a* £ A and b* £ B, 
where A, B are nonlinear sub- manifolds of the signal space M . We are given access to noisy linear 
measurements of x*: 

z = *(a* +b*) + e, (1) 

where 3? £ jg>Afx7V j g measurement matrix. Our objective is to recover the pair of signals 
(a*,b*), and thus also x*, from z. At the outset, numerous obstacles arise while trying to solve 
(1), some of which appear to be insurmountable: 
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1. (Identifiability I) Consider the simplest case, where the measurements are noiseless and the 
measurement operator is the identity, i.e., we observe x 6 such that 

x = a*+b*, (2) 

where a* £ A, b* £ B. This expression for x contains 2N unknowns but only N observa- 
tions and hence is fundamentally ill-posed. Unless we make additional assumptions on the 
geometric structure of the component manifolds A and B, a unique decomposition of x into 
its constituent signals (a*,b*) may not exist. 

2. (Identifiability II) To complicate matters, in more general situations the linear operator <1> 
in (1) might have fewer rows that columns, so that M < N. Thus, $ possesses a nontrivial 
nullspace. Indeed, we are particularly interested in cases where M <C N, in which case the 
nullspace of is extremely large relative to the ambient space. This further obscures the 
issue of identifiability of the ordered pair (a*,b*), given the available observations z. 

3. (Nonconvexity) Even if the above two identifiability issues were resolved, the manifolds A, B 
might be extremely nonconvex, or even non-differentiable. Thus, classical numerical methods, 
such as Newton's method or steepest descent, cannot be successfully applied; neither can the 
litany of convex optimization methods that have been specially designed for linear inverse 
problems with certain types of signal priors [1,6]. 

In this paper, we propose a simple method to recover the component signals (a*,b*) from z in (1). 
We dub our method Successive Projections onto Incoherent manifolds (SPIN) (see Algorithm 1) 
Despite the highly nonconvex nature of the problem and the possibility of underdetermined mea- 
surements, SPIN provably recovers the signal components (a*,b*). For this to hold true, we will 
require that (i) the signal manifolds A,B are incoherent in the sense that the secants of A are 
almost orthogonal to the secants of B; and (ii) the measurement operator 3> satisfies a certain 
restricted isometry property (RIP) on the secants of the direct sum manifold C = A® B. We will 
formally define these conditions in Section 2. We prove the following theoretical statement below 
in Section 3. 

Theorem 1 (Signal recovery) Let A,B be incoherent manifolds in~R N . Let & be a measurement 
matrix that satisfies the RIP on the direct sum C = A®B . Suppose we observe linear measurements 
z = <&(a* + b*), where a* € A and b* € B. Then, given any precision parameter v > 0, there 
exists an iterative algorithm that outputs a sequence of iterates (a^, b^) 6 A x B, k = 1, 2, . . . such 
that max{||afc — a*|| , ||b& — b*||} < Cv for all k > T. Here C and T are global positive constants. 

Our proposed algorithm (SPIN) is iterative in nature. Each iteration consists of three steps: 
computation of the gradient of the error function ^(a, b) = \ ||z — 3>(a + b)|| 2 , forming signal 
proxies for a and b, and orthogonally projecting the proxies onto the manifolds A and B. The 
projection operators onto the component manifolds play a crucial role in algorithm stability and 
performance; some manifolds admit stable, efficient projection operators while others do not. We 
discuss this in detail in Section 3. Additionally, we demonstrate that our proposed algorithm is 
stable to measurement noise (the quantity e in (1)), as well as numerical inaccuracies (such as 
finite precision arithmetic). Interestingly, SPIN exhibits a convergence rate comparable to several 
state-of-the-art algorithms [9, 10], despite the nonlinear and nonconvex nature of the reconstruction 
problem. 
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The core essence of our proposed approach has been extensively studied in a number of dif- 
ferent contexts. Methods such as Projected Landweber iterations [11], iterative hard thresholding 
(IHT) [9], and singular value projection (SVP) [10] are all instances of the same basic framework. 
SPIN subsumes and generalizes these methods. In particular, it is an iterative projected gradient 
method with the same basic approach as two recent signal recovery algorithms — Gradient Descent 
with Sparsification (GraDeS) [12], and Manifold Iterative Pursuit (MIP) [13]. We generalize these 
approaches to situations where the signal of interest is a linear mixture of signals arising from a 
pair of nonlinear manifolds. Due to this particular structure of our setting, SPIN consists of two 
projection steps (instead of one), and the analysis is more involved (see Section 4). We also explore 
the interplay between the geometric structure of the component manifolds, the linear measurement 
operator, and the stability of the recovery algorithm. 

SPIN exhibits a strong geometric convergence rate comparable to many state-of-the-art first- 
order methods. We duly note that for the case of certain special manifolds, sophisticated higher- 
order methods with stronger stability guarantees have been proposed (e.g., approximate message 
passing (AMP) [14] for sparse signal recovery and augmented Lagrangian multiplier (ALM) meth- 
ods for low-rank matrix recovery [15]). However, an appealing feature of SPIN is its conceptual 
simplicity plus its ability to generalize to mixtures of arbitrary nonlinear manifolds, provided these 
manifolds satisfy certain geometric assumptions. 

2 Geometric Assumptions 

The analysis and proof of accuracy of SPIN (Algorithm 1) requires three core ingredients: (i) a 
geometric notion of manifold incoherence that crystallizes the approximate orthogonality between 
secants of submanifolds of R. N ; (ii) a restricted isometry condition satisfied by the measurement 
operator $ over the secants of a submanifold; and (iii) the availability of projection operators that 
compute the orthogonal projection of any point x 6 R N onto a submanifold of R N . 

2.1 Manifold incoherence 

In linear inverse problems such as sparse signal approximation and compressive sensing, the as- 
sumption of incoherence between linear subspaces, bases, or dictionary elements is common. We 
introduce a nonlinear generalization of this concept. 

Definition 1 Given a manifold A C l w , a normalized secant, or simply, a secant, u G 1^ of A 
is a unit vector such that 

a a , . , 

u = — , a, a 6 A, a ^ a . 

||a — a'|| 

The secant manifold S(A) is the family of all unit vectors u generated by pairs a, a' belonging to 
A. 

Definition 2 Suppose A,B are submanifolds are ~M. N . Let 

sup |(u,u')|=e, (3) 

u&S(A), u'eS(B) 

where S(A),S(B) are the secant manifolds of A,B respectively. Then, A and B are called e- 
incoherent manifolds. 
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By definition, the quantity e is always positive; further, e < 1, due to the Cauchy-Schwartz in- 
equality. We prove that any signal x belonging to the direct sum A®B can be uniquely decomposed 
into its constituent signals when the upper bound on e holds with strict inequality. 

Lemma 1 (Uniqueness) Suppose that A,B are e-incoherent with < e < 1. Consider x = 
a + b = a' + b', where a, a' £ A and b, b' £ B. Then, a = a', b = b'. 

Proof. It is clear that ||a + b — (a' + b')|| 2 = 0, i.e., 

||a - a'f + ||b - b'|| 2 = -2(a - a',b - b') < 2 |(a - a', b - b') | . 

However, due to the manifold incoherence assumption, the (unnormalized) secants a — a', b — b' 
obey the relation: 

|<a-a',b-b')| <e||a-a'|| ||b-b'|| < |e(||a - a'|| 2 + ||b - b'f), (4) 

where the last inequality follows from the relation between arithmetic and geometric means (hence- 
forth referred to as the AM-GM inequality). Therefore, we have that 



'|| 2 + ||b-b'|| 2 <e(||a-a'|| 2 + ||b-b'|| 2 ) , 



||a — a 

for e < 1, which is impossible unless a = a', b = b'. □ 
We can also prove the following relation between secants and direct sums of signals lying on 
incoherent manifolds. 

Lemma 2 Suppose that A,B are e-incoherent with < e < 1. Consider xi = ai + bi,X2 = a2 + b2 
where ai, a2 £ A and bi, b2 £ B. Then 



|(ai - a 2 ,bi - b 2 )| < 2 {\-e) ^ ~ X2 ^ 



Proof. From (4), we have 



|(ai-a 2 ,bi-b 2 )| < |(||ai - a 2 || 2 + ||bi - b 2 || 2 ) 

e 2 

= ^ H ai + bl ~ &2 ~ b2 'l ~ e ( a l ~ a 2> b l - b 2) 

< -||xi-x 2 || +e|(ai - a 2 ,bi - b 2 )| . 
Rearranging terms, we obtain the desired result. □ 



2.2 Restricted isometry 

Next, we address the situation where the measurement operator $ G R MxN contains a nontrivial 
nullspace, i.e., when M < N. We will require that $ satisfies a restricted isometry criterion on the 
secants of the direct sum manifold C = A © B. 

Definition 3 Let C be a submanifold ofM. N . Then, the matrix $ 6 M MxAr satisfies the restricted 
isometry property (RIP) on C with constant 5 £ [0,1), if for every normalized secant u belonging 
to the secant manifold S(C), we have that 

1-S< ||$u|| 2 < l + <5. (5) 
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The notion of restricted isometry (and its generalizations) is an important component in the 
analysis of many algorithms in sparse approximation, compressive sensing, and low-rank matrix 
recovery [4, 6]. While the RIP has traditionally been studied in the context of sparse signal models, 
(5) generalizes of this notion to arbitrary nonlinear manifolds. The restricted isometry condition is 
of particular interest when the range space of the matrix <1> is low-dimensional. A key result [16] 
states that, under certain upper bounds on the curvature of the manifold C, there exist probabilistic 
constructions of matrices <1? that satisfy the RIP on C such that the number of rows of is 
proportional to the intrinsic dimension of C, rather than the ambient dimension N of the signal 
space. We will discuss this further in Section 5. 

2.3 Projections onto manifolds 

Given an arbitrary nonlinear manifold A 6 R^, we define the operator V_a( ) as the Euclidean 
projection projector onto A: 

Va(x) = arg min ||x' — x|| . (6) 

We observe that for arbitrary nonconvex manifolds A, the above minimization problem (6) may 
not yield a unique optimum. Technically, therefore, "P/i(x) is a set- valued operator. For ease of 
exposition, V^,{x) will henceforth refer to any arbitrarily chosen element of the set of signals that 
minimize the ^-error in (6). 

The projection operator V_a plays a crucial role in the development of our proposed signal 
recovery algorithm in Section 3. Note that in the case of general nonlinear manifolds, Va may 
be quite difficult to compute exactly. Therefore, following the lead of [13], we also define a 7- 
approximate projection operator onto A: 

x' = V\(x) =>■ x' £ A, and ||x' - x|| 2 < ||^(x) - xf + 7, (7) 

so that V\{x) yields a vector x' £ A that approximately minimizes the squared distance from x to 
A. Again, 7^(x) need not be uniquely defined for a particular x. 

Certain specific instances of nonconvex manifolds do admit efficient exact projection operators. 
For example, consider the space of all A-sparse signals of length N; this can be viewed as the 
union of the (^) canonical subspaces in R^ or, alternately, a A-dimensional submanifold of R w . 
Then, the projection of an arbitrary vector x £ R w onto this manifold is merely the best A-sparse 
approximation to x, which can be very efficiently computed via simple thresholding. We discuss 
additional examples in Section 5. 

3 The SPIN Algorithm 

We now describe an algorithm to solve the linear inverse problem (1). Our proposed algorithm, 
Successive Projections onto IN coherent manifolds (SPIN), can be viewed as a generalization of 
several first order methods for signal recovery for a variety of different models [9, 10, 13]. SPIN is 
described in pseudocode form in Algorithm 1. 

The key innovation in SPIN is that we formulate two proxy vectors for the signal components 
a/c and and project these onto the corresponding manifolds A and B. We demonstrate that 
SPIN possesses strong uniform recovery guarantees comparable to existing state-of-the-art algo- 
rithms for sparse approximation and compressive sensing, while encompassing a very broad range 
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Algorithm 1 Successive Projections onto Incoherent Manifolds (SPIN) 

Inputs: Observation matrix <E>, measurements z, projection operators V^(-),Vb(-), 

number of iterations T, step size 77 
Outputs: Estimated signal components aei,bG6 
Initialize: ao = 0, bo = 0, r = z, k = 



while k < T do 

g fc <- $ T r 

a fe <- a fc + ?7g fc , b fc ^ b fc + r/gfc 
afe+i <- £U(3fc), b fe+ i <- PB(bfe) 
r <- z - $(a fc+ i + b fe+ i) 

k <r- k + 1 

end while 

return (a, b) <— (ay, by) 



{form gradient} 
{form signal proxies} 
{apply projection operators} 
{update residual} 



of nonlinear signal models. The following theoretical result describes the performance of SPIN for 
signal recovery. 

Theorem 2 (Main result) Suppose A,B are e-incoherent manifolds in R N . Let <& be a measure- 
ment matrix with restricted isometry constant 5 over the direct sum manifold C = A © B. Suppose 
we observe noisy linear measurements z = $(a* + b*) + e, where a* £ A and h* £ B. If 

1 — 9e 

< 5 < - — — , (8) 
3 + 5e' V ' 

then SPIN (Algorithm 1) with step size 77 = 1/(1 + 5) with exact projections Va-i'Pb outputs &t G A 
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and hr £ B, such that ||z — $(a^ + br)|| < f3 ||e|| + v in no more than T = \ i og ^/ a ^ log 2it1 
iterations for any v > 0. Here, a,f3 are moderately sized positive constants. 

For the special case when there is no measurement noise (i.e., e = 0), Theorem 2 states that 
after a finite number of iterations, SPIN outputs signal component estimates (a, b) such that 

z — $(a + b) < v for any desired precision parameter v. From the restricted isometry assumption 
on $ and Lemma 1, we immediately obtain Theorem 1. Since we can set v to an arbitarily 
small value, we have that the SPIN estimate (a, b) converges to the true signal pair (a*,b*). 
Exact convergence of the algorithm might potentially take a very large number of iterations, but 
convergence to any desired positive precision constant (3 takes only a finite number of iterations. 
For the rest of the paper, we will informally denote signal "recovery" to imply convergence to a 
sufficiently fine precision. 

SPIN assumes the availability of the exact projection operators Vj\,,Vb- In certain cases, it 
might be feasible to numerically compute only 7-approximate projections, as in (7). In this case, 
the bound on the norm of the error z — $(a^ + br) is only guaranteed to be upper bounded by 
a positive multiple of the approximation parameter 7. The following theoretical guarantee (with a 
near-identical proof mechanism as Theorem 2) captures this behavior. 

Theorem 3 (Approximate projections) Under the same suppositions of Theorem 2, SPIN (Al- 
gorithm 1) with ^-approximate projections and step size 77 = 1/(1 + 5) outputs a? £ A and by € B 
such that ||z — <&(ax + br)|| 2 < (3 \\e\\ 2 + j^^f+u , in no more thanT = [~ log ^/ a ) log ^j-] iterations. 
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We note some implications of Theorem 2. First, suppose that <3? is the identity operator, i.e., 
we have full measurements of the signal x* = a* + b*. Then, 5 = and the lower bound on 
the restricted isometry constant holds with equality. However, we still require that e < 1/9 for 
guaranteed recovery using SPIN. We will discuss this condition further in Section 5. 

Second, suppose that the one of the components manifolds is the trivial (zero) manifold; then, 
we have that e = 0. In this case, SPIN reduces to the Manifold Iterative Pursuit (MIP) algorithm for 
recovering signals from a single manifold [13]. Moreover, the condition on 5 reduces to < 5 < 1/3, 
which exactly matches the condition required for guaranteed recovery using MIP. 

Lastly, the condition (8) in Theorem 2 automatically implies that e < 1/9. This represents a 
mild tightening of the condition on e required for a unique decomposition (Lemma 1), even with 
full measurements (i.e., when <1? is the identity operator or, more generally, when (5 = 0). 



4 Analysis 

The analysis of SPIN is based on the proof technique developed by [12] used in analyzing the 
iterative hard thresholding (IHT) algorithm for sparse recovery and further extended in [10,13]. 
For a given set of measurements z obeying (1), define the error function ijj : A x B — > R as 

^(a,b) = i||z-*(a + b)|| 2 . 

It is clear that ip(a*,h*) = \ ||e|| 2 . The following lemma bounds the error of the estimated signals 
output by SPIN at the (k + l)-th iteration in terms of the error incurred at the k-th iteration, and 
the norm of the measurement error. 

Lemma 3 Define (afc,bfc) as the intermediate estimates obtained by SPIN at the k-th iteration. 
Let 5, e be as defined in Theorem 2. Then, 

V>(afc+i, bfc+i) < aip(a k , bfc) + C ||e|| 2 , (9) 

where 

_2^_ _i_ <=;!+i_^ l \_A l±§. 6 

1-8 1 — 8 l-e p 2 1 — S l-e 

1_9i+i_^ ' 1 n l+8 e ■ 

1 z 1-6 l-e 1 * 1-6 l-e 

Proof. Fix a current estimate of the signal components (a&, b&) at iteration k. Define x& = a^ + b^. 
Then, for any other pair of signals (a, b) € A x B, we have 



V(a,b)-V(a fc ,b fc ) = i(||z-*(a + b)|| 2 -||z-*(a fc + b fc )|| 2 ) 



= I||z-*(a + b)|| 2 -i||z-*(a + b)|| 2 + (*(a fc + b fc )-*(a + b), z) 

1 2 

= -||*x-*x fc || + (z - *x fc ,*x fc - *x), 

where x = a + b, and the last equation is obtained by completing the square. Since $ is a linear 
operator, we can take the adjoint within the inner product to obtain 

1 



V>(a,b) - V(a fc ,b fc ) = ^ ||*x - $x fe || 2 + ($ T (z - *x fc ), x fc - x) 

1 



< ^(l + <5)||x-x fe || 2 + ($ T (z-$x fc ),x fc -x). (10) 
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The last inequality occurs due to the RIP of * applied to the secant vector x — Xfc G <S(C). To 
the right hand side of (10), we further add and subtract 2 (i+s) ||^ T ( Z ~~ ^ x fc)|| 2 to complete the 
square: 



V>(a,b) - V(a fc ,b fe ) < ^(1 + 6) 



1 Tl 
l + o 



' ]* T (z-*x fe ) 



2(1 + 5) 

Define gfc = yq^3? T ( z — &(a k + bfc)) as the gradient of the error function at (a&, bfc). Then, 

^(a,b)-^(a*,b fc ) < + (l|a + b-(a fc + b fc + g fc )|| 2 - ||g fe || 2 ) . (11) 

Next, define the function ( on A x B as £(a, b) = ||a + b — (afc + bfc + gfc)|| 2 . Then, we have 

C(a fe+ i,b fe+ i) = ||a fc+1 - (a fc + g fc ) + b fe+ i - (b fc + g fe ) + g fc || 2 

= ||a fc+ i - (a fe + g fc )|| 2 + ||b fe+ i - (b fe + g fc )|| 2 + ||g fc || 2 

+2(a fc+1 - (a fe + g fc ),b fc+1 - (b fe + g fc )) + 2(g fc ,a fc+1 + b fc+1 - (a fc + h k + 2g fc )). 

But, by definition, a fc+ i = Vj^k + gfc), and hence ||a fe+ i - (a k + gfc) II < ||a - (a* + gfc)|| for any 
a G A. An analogous relation can be formed between bfc +1 and b*. Hence, we have 

||a fe+ i - (a^ + gfc)|| < ||a* - (a fc + g fc )|| , and 

||bfc+i-(bfc + gfc)|| < ||b*-(b fc + gfc)||. 

Substituting for (afc+i,bfc+i), we obtain 

C(afc+i,bfc+i) < ||a*-(a fc + g fc )|| 2 + ||b*-(b fc + gfc)|| 2 + ||gfc|| 2 

+2(a fc+ i - (a fe + gfc),b fc+ i - (b fc + gfc)) + 2(g fc ,afc + i + b fc+ i - (a fe + b fc + 2g fc )) 
= ||a* - (afc + g fc )|| 2 + ||b* - (bfc + gfc)|| 2 + ||gfc|| 2 

+2(a* - (afc + g fe ), b* - (bfc + gfc)) + 2(gfc, a* + b* - (a fc + b k + 2g fe )) 
+2(afc +1 - (afc + gfc),bfc +1 - (b fc + gfc)) - 2(a* - (a* + g fc ),b*. +1 - (b fe + g fe )) 
+2(gfc,a fe+ i + bfc +1 - (a* + b*)). 

Completing the squares, we have: 

C( a fc+i,bfc + i) < ||a* + b* - (a fc + b fc + g fc )|| 2 + 2(afc +i - afc,b fc+ i - b fe ) - 2(a* - a fc ,b* - b fe ) 

+2(gfc, -afc+i + a fc - b fe+ i + b fc + a* - a fc + b* - b fe + a k+1 + b fc+1 - (a* + b*)). 

The last term on the right hand side equals zero, and so we obtain 

C(a*+i,b fc+ i) < C(a*,b*) + 2(a fe+ i - afc,b fc+ i - b fc ) - 2(a* - a fe ,b* - bfc). 

Combining this inequality with (11), we obtain the series of inequalities 

ip(a k+1 ,b k+1 ) - V(a fc ,bfc) < -(1 + 6) (c(afc+i> h k+1 ) - ||gfc|| 2 ) 

Ti 



< I(l + 5)( C (a*,b*)-||gfc|| 2 ) 



T 2 



+ (1 + 5) ((afc+i - afc, bfc + i - bfc) - (a* - a k ,h* - h k )) (12) 
= Ti + T 2 . 



8 



We can further bound the right hand side of (12) as follows. First, we expand Ti to obtain 

Ti = ^(l + <5)(||a* + b*-(a fc + b fc + g fc )|| 2 -||g fc || 2 ) 

= 1(1 + 5) (||a* + b* - (a* + b fc )|| 2 + 2<g fc , a* + b* - (a* + b fe )>) 
= -(1 + S) ||x* - x fc || 2 + ($ T (z - *x fc ),x fc - x*). 
Again, x* — x& is a secant on the direct sum manifold C. By the RIP property of we have 

Ti < i^||*x*-*x fe || 2 + ($ T (z-*x A; ),x fc -x*) 

= X - ||*x* - *x fc || 2 + (* T (z - *x fc ),x fc - x*) + ||z - $x fc || 2 

2(5 

= V(a*,b*)-^(a fc ,b fc ) + ^— -^(a fc ,b fc ), (13) 

where the last expression follows easily from (10) and the definition of ip. 
The term T2 can be bounded using Lemma 2 as follows. We have 

(a* - a fc ,b* - b fe ) < - ||x* - x fc || 2 , and 

(afe+i - a k ,b k+ i - bfe) < — - ^ ||x fc+ i - x fc || 2 

^ (I^) (l|x fe+ i-x*|| 2 + ||x fc -x*|| 2 ), (14) 

where (14) follows from the the triangle inequality and the AM-GM inequality. Further, T2 can be 
simplified as 

T 2 < (1 + 5)^ Q||x*-x fe || 2 + ||x*-xfe +1 || 2 



But, 



|*x*-*x fc || 2 = ||$x* - *x fe + e - e|| 2 



< 2 || *x* - *Xfe + e|| 2 + 2 ||e|| 2 = 4^(a fe , b fc ) + 4^(a*, b*), 
via the same technique used to obtain (14). Similarly, 

|| *x* - *x fc+1 || 2 < 4^(a*+i,b fc+ i) + 4^(a*,b*). 

Hence, we obtain 

T 2 < ] + S . n 26 (3^>(ajfe, bfe) + ^(ajfe+i, b fc+ i) + 4^(eu, b>)) ■ (15) 
1 — 1 — e 
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Combining (12), (13), and (15), we obtain 



V>(a fe+1 ,b fe+1 ) < ^(a*,b*) + T ^V(a fe ,b fc ) 

+137 ( 3 ^ afc ' bfc ) + b fc+ i) + 4^(a*, b*)) . 

Rearranging, we obtain Lemma 3. □ 
Proof of Theorem 2. Equation 9 describes a linear recurrence relation for the sequence of positive 
real numbers Y>(afc, b^), k = 0, 1,2, . . . with leading coefficient a < 1. By choice of initialization, 

V>(ao,bo) = Therefore, for all fc 6 N, we have the relation 

^(afc, b fe ) < a fe V(a , b ) + ||e|| 2 . 

1 — a 

II 1 1 2 

Choosing /3 = yz^, and k > T such that T = r iog(i/a) log ^rl j the result follows. □ 



5 Applications 

The two-manifold signal model described in this paper is applicable to a wide variety of problems 
that have attracted considerable interest in the literature over the last several years. We discuss a 
few representative instances and show how SPIN can be utilized for efficient signal recovery in each 
of these instances. We also present several numerical experiments that indicate the kind of gains 
that SPIN can offer in practice. 



5.1 Sparse representations in pairs of bases 

We revisit the classical problem of decomposing signals in an overcomplete dictionary that is the 
union of a pair of orthonormal bases and show how SPIN can be used to efficiently solve this 
problem. Let A, A' be orthonormal bases of M N . Let A be the set of all ^-sparse signals in M. N 
in the basis expansion of A, and let B be the set of all K2- sparse signals in the basis expansion of 
A'. Then, A and B can be viewed as K\- and ^-dimensional submanifolds of W N , respectively. 
Consider a signal x* = a* + b*, where a*ei,b*eB so that 

i=i i=i 

The problem is to recover (a*, b*) given x*. This problem has been studied in many different forms 
in the literature, and several algorithms have been proposed in order to solve it efficiently [1, 3, 17]. 
See [18] for an in-depth study of the various state-of-the-art methods. All these methods assume 
a certain notion of incoherence between the two bases, most commonly referred to as the mutual 
coherence /j,, which is defined as 

//(A, A') = max |(Aj, Aj)| . (16) 

It is clear that fj, < 1. It can be shown that the mutual coherence also obeys the lower bound 
H > \/y/~N for any basis of M. N . This lower bound is in fact tight; equality is achieved, for example, 
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when A is the canonical basis in R^ and A' is the basis defining the Walsh-Hadamard transform 
or the Discrete Fourier Basis [3, 17]. 

We establish the following simple relation between ^ and the manifold incoherence between A 
and B. 

Lemma 4 Let A be the set of all K\-sparse signals in K, and B be the set of all K^-sparse signals 
in A'. Let e denote the manifold incoherence between A and B. Then, 

e<n{K,M){K l +K 2 ). 

Proof. The secant manifold S(A) is equivalent to the set of signals in W N that are 2K\ -sparse in 
A; similarly, S(B) is equivalent to the set of 2i^2-sparse signals in B. Therefore, if one considers 
unit norm vectors u € S(A), u' G S(B), we obtain 



u,u' = 



/2K! 2K 2 



2Kt 2K 2 



2Ki 2K 2 

i j 



/2K X \ l2K 2 



(17) 



Kl = l 



of [17], we have 
l,...,2Ki and 



where the last relation follows from the triangle inequality. From Theorem 1 
that the right hand side is maximized when a, = ' ^ = \/2K 2 ^ or a ^ * 
j = 1, . . . , 2K 2 . The result of Lemma 4 follows from inserting these values into (17) and applying 
the AM-GM inequality. □ 
We show how SPIN can be used to solve the linear inverse problem of recovering (a* , b* ) 
from x*. The restricted isometry assumption is not relevant in this case, since we assume that 
we have full measurements of the signal; therefore 5 = 0. An upper bound for the manifold 
incoherence parameter e is specified in Lemma 4. The (exact) projection operators Va,Vb can be 
easily implemented; we simply perform a coefficient expansion in the corresponding orthornormal 
basis and retain the coefficients of largest magnitude. Mixing these ingredients together, we can 
guarantee that, given any signal x*, SPIN will return the true components a*,b*. This guarantee 
is summarized in the following result. 

Corollary 1 (SPIN for pairs of bases) Let x* = a* + b* ; where a* is K\-sparse in A and b* 
is K-i-sparse in A'. Let [i denote the mutual coherence between A and A'. Then, SPIN exactly 
recovers (a*,h*) from x* provided 

(18) 

Proof. If (18) holds, then from Lemma 4 we know that the manifold incoherence e between A and 
B is smaller than 1/9. But this is exactly the condition required for guaranteed convergence of 
SPIN to the true signal components (a*,b*). □ 
SPIN thus offers a conceptually simple method to separate mixtures of signals that are sparse 
in incoherent bases. The best known approach for this problem is an t\ -minimization formulation 
that features a similar guarantee that is known to be tight [17, 19]: 



K 1 +K 2 < 



v/2-0.5 



/i 
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Therefore, SPIN yields a recovery guarantee that is equivalent to the best possible method, modulo 
a small constant. It is likely that the constant 1/9 in Corollary 1 can be sharpened, but we will 
not pursue this direction in this paper. 



5.2 Articulation manifolds 

Articulation manifolds provide a powerful, flexible conceptual tool for modeling signals and image 
ensembles in a number of applications [20, 21]. Consider an ensemble of signals At C R^ that are 
generated by varying K parameters 9 £ 0, G C R . Then, we say that the signals trace out a 
nonlinear i^-dimensional articulation manifold in WL N , where 9 is called the articulation parameter 
vector. Examples of articulation manifolds include: acoustic chirp signals of varying frequencies 
(where 9 represents the chirp rate) ; images of a white disk translating on a black background (where 
9 represents the planar location of the disk center); and images of a solid object with variable pose 
(where 9 represents the 6D pose parameters, three corresponding to spatial location and three 
corresponding to orientation). 

We consider the class of compact, smooth, K-dimensional articulation manifolds Ai C M. N . For 
such manifold classes, it is possible to construct linear measurement operators 3> that preserve the 
pairwise secant geometry of Ai. Specifically, it has been shown [16] that there exist randomized 
constructions of measurement operators 6 R MxiV that satisfy the RIP on the secants of Ai with 
constant 8, and with probability at least p, provided 



for some constant C that depends only on the smoothness and volume of the manifold Ai. There- 
fore, the dimension of the range space of 3> is proportional to the number of degrees of freedom K, 
but is only logarithmic in the ambient dimension N. Moreover, given such a measurement matrix 
$ with isometry constant 5 < 1/3 and a projection operator Vm(') onto Ai, any signal x € Ai 
can be reconstructed from its compressive measurements y = 3>x using Manifold Iterative Pursuit 



We generalize this setting to the case where the unknown signal of interest arises as a mixture 
of signals from two manifolds A and B. For instance, suppose we are interested in the space of 
images, where A and B comprise of translations of fixed template images /(t) and g(t), where t 
denotes the 2D domain over which the image is defined. Then, the signal of interest is an image of 
the form 



where 9\ and #2 denote the unknown translation parameters. The problem is to recover (a*,b*), 
or equivalently (9±, O2), given compressive measurements z = <l>(a* + b*). 

We demonstrate that SPIN offers an easy, efficient technique to recover the component images. 
This example also demonstrates that SPIN is robust to practical considerations such as noise. 
Figure 1 displays the results of SPIN recovery of a 64 x 64 image from very limited measurements. 
The unknown image consists of the linear sum of arbitrary translations of template images /(t) and 
g(t), that are smoothed binary images on a black background of a white disk and a white square, 
respectively. Further, the image has been contaminated with significant Gaussian noise (SNR = 
14dB) prior to measurement (Fig. 1(a)). From Figs. 1(b) and 1(c), we observe that SPIN is able to 
perfectly recover the original component signals from merely M = 50 random linear measurements. 




(MIP) [13]. 



x* =a*+b* = /(t + 0i) + 5 (t + 2 ) 
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(a) Original image 



(b) Disk component 



(c) Square component 



Figure 1: SPIN recovery of a noisy 64 x 64 image from compressive measurements. The clean image consists 
of the linear superposition of a disk and a square of fixed pre-specihed sizes, but the locations of the centers 
of the disk and the square are unknown. Additive Gaussian noise (SNR = 14dB) is added to the image prior 
to measurement. Signal length N — 64 x 64 = 4096, number of compressive measurements M = 50. (a) 
Original noisy image, (b) Reconstructed disk, (c) Reconstructed square. SPIN perfectly reconstructs both 
components from just M/N = 1.2% measurements. 

For guaranteed SPIN convergence, we require that the manifolds A, B are incoherent. Infor- 
mally, the condition of incoherence on the secants of A and B is always valid when the template 
images f(t),g(t) are "sufficiently" distinct. This intuition is made precise using the bound in (3). 
More generally, we can state the following theoretical guarantee for SPIN performance in the case 
of general higher-dimensional manifolds. 

Corollary 2 (SPIN for pairs of manifolds) Let the entries of $ G l MxJV be chosen from a 
standard Gaussian probability distribution. Let A and B be e-incoherent compact submanifolds of 
~R N of dimensions K and K' respectively. Let z = $(a* + b*), where a* 6 A and b* 6 B. Then, 
with probability greater than 1 — p, SPIN exactly recovers (a*, b*) from z, provided 



Here, Ca,Cb are constants that depend only on certain intrinsic geometric parameters (such as 
volume) of A,B respectively. 

Proof. It is easy to show that the matrix $ satisfies the RIP on the secants of the direct sum 
C = A® B if (19) is true. This can be proved using a simple union bound argument in conjunction 
with the techniques used in Section 3.2.2 in [16]. Once the RIP of <& is established, SPIN recovery 
follows from Theorem 2. □ 
An important consideration in SPIN is the tractable computation of the projection Vm( x ) 
given any x 6 M N . For example, in the numerical example in Fig. 1, the operator V^{x) onto the 
manifold A consists of running a matched filter between the template f(t) and the input signal x 
and returning f(t + 9), where the parameter value 9 corresponds to the 2D location of the maximum 
of the matched filter response. This is very efficiently carried out in O (N log N) operations using 
the Fast Fourier Transform (FFT). However, for more complex articulation manifolds, the exact 
projection operation might not be tractable, whereas only approximate numerical projections can 
be efficiently computed. In this case also, SPIN can recover the signal components (a*,b*), but 
with weaker convergence guarantees (Theorem 3). 



M > O ((K\og(C A N) + K'\og{C B N)) log^ 1 )) . 



(19) 



13 



(a) Original signal 



(b) Manifold component 



(c) Noise component 



Figure 2: SPIN recovery of a shifted Gaussian pulse from compressive measurements. The shift parameter of 
the pulse is unknown, and the signal is corrupted with K' -sparse, impulsive noise of unknown amplitudes and 
locations. N — 10000, 7^' = 10, M — 150. (a) Original signal, (h) Reconstructed Gaussian pulse (Recovery 
SNR — 80.09 dB). (c) Estimated noise component. SPIN perfectly reconstructs both components from just 
M/N = 1.5% measurements. 

5.3 Signals in impulsive noise 

In some situations, the signal of interest x might be corrupted with impulsive noise (or shot noise) 
prior to signal acquisition via linear measurements. For example, consider Fig. 2(a), where the 
Gaussian pulse is the signal of interest, and the spikes indicate the undesirable noise. In this case, 
the linear observations are more accurately modeled as: 

z = <&(x + n), such that x € M, 

and n is a i^'-sparse signal in the canonical basis. Therefore, SPIN can be used to recover x from 
z, provided that the manifold M. is incoherent with the set of sparse signals S^/ and 3? satisfies 
the RIP on the direct sum Ai + X^'. 

Figure 2 displays the results of a numerical experiment that illustrates the utility of SPIN in 
this setting. We consider a manifold of signals of length N = 10000 that consist of shifts of a 
Gaussian pulse of fixed width <?o £ R . The unknown signal x is an element of this manifold A4, 
and is corrupted by K' = 10 spikes of unknown magnitudes and locations. This degraded signal is 
sampled using M = 150 random linear measurements to obtain an observation vector z. 

We apply SPIN to recover x from z. The projection operator Vj^i') consists of a matched filter 
with the template pulse go, while the projection operator Vs K ,(-) simply returns the best K'-tevm 
approximation in the canonical basis. Assuming that we have knowledge of the number of nonzeros 
in the noise vector n, we can use SPIN to reconstruct both x and n. We observe from Fig. 2(b) 
that SPIN recovers the true signal x with near-perfect accuracy. Further, this recovery is possible 
with only a small number M = 150 linear measurements of x, which constitutes but a fraction of 
the ambient dimension of the signal space. 

Figure 3 plots the number of measurements M vs. the signal reconstruction error (normalized 
relative to the signal energy, and plotted in dB). We observe that, by increasing M, SPIN can 
tolerate an increased number K' of nuisance spikes. Further, by Corollary 2, we observe that this 
relationship between M and K' is in fact linear. This result can be extended to any situation where 
the signals of interest obey a "hybrid" model that is a mixture of a nonlinear manifold and the set 
of sparse signals. 



14 




Figure 3: Monte Carlo simulation of SPIN signal recovery in impulsive noise , averaged over 100 trials. In 
each trial, the measured signal is the sum of a randomly shifted Gaussian pulse and a random K 1 -sparse 
signal. SPIN can tolerate a higher number K' nuisance impulses by increasing the number of measurements 
M. By Corollary 2 this dependence of M on K' can be shown to be linear. 

6 Discussion 

We have proposed and rigorously analyzed an algorithm, which we dub Successive Projections 
onto INcoherent Manifolds (SPIN), for the recovery of a pair of signals given a small number of 
measurements of their linear sum. For SPIN to guarantee signal recovery, we require two main 
geometric criteria to hold: (i) the component signals should arise from two disjoint manifolds that 
are in a specific sense incoherent, and (ii) the linear measurement operator should satisfy a restricted 
isometry criterion on the secants of the direct sum of the two manifolds. The computational 
efficiency of SPIN is determined by the tractability of the projection operators onto either component 
manifold. We have presented indicative numerical experiments demonstrating the utility of SPIN, 
but defer a thorough experimental study of SPIN to future work. 

Practical considerations. SPIN is an iterative gradient projection algorithm and requires as 
input parameters the number of iterations T and the gradient step size rj. The iteration count T 
can be chosen using one of many commonly-used stopping criteria. For example, convergence can 
be declared if the norm of the error b&) at the {k + l)-th time step does not differ significantly 
from the error at k-th time step. The choice of optimal step size rj is more delicate. Theorem 2 
relates the step size to the restricted isometry constant d of but this constant is not easy to 
calculate. In our preliminary findings, a step size in the range 0.5 < rj < 0.7 consistently gave good 
results. See [22] for a discussion on the choice of step size for hard thresholding methods. 

In practical scenarios, the signal of interest rarely belongs exactly to a low-dimensional subman- 
ifold Ai of the ambient space, but only is well-approximated by A4. Interestingly, in such situations 
the effect of this mismatch can be studied using the concept of 7-approximate projections (7). The- 
orem 3 rigorously demonstrates that SPIN is robust to such approximations. Further, our main 
result (Theorem 2) indicates that SPIN is stable with respect to inaccurate measurements, owing to 
the fact that the reconstruction error is bounded by a constant times the norm of the measurement 
noise vector e. 

More than two manifolds. For clarity and brevity, we have focused our attention on signals 
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belonging to the direct sum of two signal manifolds. However, SPIN (and its accompanying proof 
mechanism) can be conceptually extended to sums of any Q manifolds. In such a scenario, the con- 
ditions of convergence of SPIN would require that the component manifolds are Q-wise incoherent, 
and the measurement operator $ satisfies a restricted isometry on the Q-wise direct sum of the 
component manifolds. 

Connections to matrix recovery. An intriguing open question is whether SPIN (or a similar 
first-order projected gradient algorithm) is applicable to situations where either of the component 
manifolds is the set of low-rank matrices. The problem of reconstructing, from affine measurements, 
matrices that are a sum of low-rank and sparse matrices has attracted significant attention in the 
recent literature [7,8,23]. The key stumbling block is that the manifold of low-rank matrices is 
not incoherent with the manifold of sparse matrices; indeed, the two manifolds share a nontrivial 
intersection (i.e., there exist low rank matrices that are also sparse, and vice versa). Phenomena 
such as these make the analysis of SPIN (or similar algorithms) quite challenging, and it may be 
that higher-order techniques for signal reconstruction will be needed. 
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